!-------------------------------------------------!
!               MODULE OUTCOMES			  !
!-------------------------------------------------!
! This module calculates all the properties       !
! of the sistem	                		  !
!-------------------------------------------------!

  module outcomes
  
    use constants
    use control
    use gcontacts
    use gdevice
    use matrixmod 
    use integration
    use hamiltonian
 
    implicit none
    
    private 
    
!-------------------------------------------------!
! Variables declaration	          		  !					
!-------------------------------------------------!
!-------------------------------------------------!
! Public subroutines	 	        	  !					
!-------------------------------------------------!
    
    public :: DOS_contact, DOS_device
    public :: current, IvsV_curve, transmision_function
    public :: IVSinglePoint, RhoSinGlePoint, GSinglePoint, current_rho, IvsL_curve
    public :: eigenvals	

  contains

!-------------------------------------------------!
! 	DOS CONTACT		 		  !					
!-------------------------------------------------!

   subroutine DOS_contact()
   
   integer :: i
   complex(dp) :: z
      
   call cpu_time(t_1)  

   open(1,file='DOS_contacts.dat')  
     
    do i=1,NEP_ContactsDOS
      z=cmplx(LEV_ContactsDOS + i*ES_ContactsDOS, RAP_constant)
      write(1,*)real(z),(-1.0_dp/pi)*aimag(G0_contact(z,0.0_dp))
    enddo
    
   close(1)
   call cpu_time(t_2)  

   write(9,*)'DOS_contact_time',t_2-t_1
	    
   end subroutine DOS_contact

!-------------------------------------------------!
! DOS DEVICE: it calculates the trace		  !
! of the G+ function device    	 		  !					
!-------------------------------------------------!

   subroutine DOS_device()
   
   integer :: i
   complex(dp) :: z, dl(N,N), dr(N,N), gp(N,N), suma(N,N), gd(N,N), Aux(N,N)
   
   call cpu_time(t_1)     
   write(*,*)'DOS device ...'  
   open(2,file='DOS_device.dat')  
     
   if(UseSMatrix_logic.EQ.'N')then
     do i=1,NEP_DeviceDOS
       z=cmplx(LEV_DeviceDOS + i*ES_DeviceDOS,RAP_constant)
       write(2,*)real(z), (-1.0_dp/pi)*aimag(trace(G_dev(z)))
     enddo      
   else
     Aux=cmplx(1.0_dp,0.0_dp)*Solap
     do i=1,NEP_DeviceDOS
       z=cmplx(LEV_DeviceDOS + i*ES_DeviceDOS,RAP_constant)
       write(2,*)real(z), (-1.0_dp/pi)*aimag(trace(G_dev(z).x.Aux))
     enddo      
   endif
    
   close(2)
 
   call cpu_time(t_2)     
   write(9,*)'DOS_device_time',t_2-t_1
	    
   end subroutine DOS_device
   
!-------------------------------------------------!
! Transmision Function 	   	 		  !					
!-------------------------------------------------!

   subroutine transmision_function()
   
   integer :: i
   complex(dp) :: z
   real(dp) :: E
   
   write(*,*)'transmision funtion ...'  
   open(1,file='transmision_function.dat')  
 
    do i=1,NEP_TransmisionF
      E=LEV_TransmisionF + i*ES_TransmisionF
      write(1,*)E,transmision(E)
    enddo
    
   close(1)
	    
   end subroutine transmision_function
     
!-------------------------------------------------!
! Electric Current with transmision function      !					
!-------------------------------------------------!

   real(dp) function current(mu1,mu2)
   
   integer :: i
   real(dp) :: E,mu1,mu2,dE
   
   dE=(mu1-mu2)/real(N_steps)
   current=0.0_dp
        
    do i=1,N_steps
      current=current + dE*transmision(mu2+dE*i)
    enddo
    
    current=(2.0_dp*mAtoefs_constant/(hbar*2.0_dp*pi))*current

   end function current

!-------------------------------------------------!
! Electric Current calculated with rho 		  !					
! This works only with S=I
!-------------------------------------------------!

   real(dp) function current_rho(mu1,mu2)
   
   integer :: i,j
   real(dp) :: E,mu1,mu2,dE
   complex(dp) :: rho_singlepoint(N,N),rp(N,N), Aux(N,N)
   
   mu_left=mu1
   mu_right=mu2
   
   current_rho=0.0_dp
    
   call initialize_integration()
   rho_SinglePoint=path_integrator()
   call finalize_integration()
   
   do i=1,N 
     do j=1,N
        current_rho=(4.0_dp*mAtoefs_constant/(hbar))*H(i,j)*aimag(rho_SinglePoint(i,j))
!        rp=rho_SinglePoint
	write(12345,*)current_rho
     enddo 
   enddo

  ! Aux=cmplx(1.0_dp,0.0_dp)*Solap
  ! write(*,*)trace(rho_SinglePoint.x.Aux)
  ! stop

   end function current_rho
   
!-------------------------------------------------!
! IvsV curve 	 	   	 		  !					
!-------------------------------------------------!

   subroutine IvsV_curve() 
      
   integer :: i
   real(dp) :: deltaV,Voltage,Voltage0,mu1,mu2
   
   Voltage0=-LVV_IvsV

   deltaV=2.0_dp*Voltage0/real(NVP_IvsV)
   
   write(*,*)'IV curve ...'
   
   open(1,file='IV_curve.dat')
   
    do i=1,NVP_IvsV
      
      write(*,*)'Vstep',i,NVP_IvsV
      Voltage=deltaV*i - Voltage0
    
      mu2=-Voltage/2.0_dp
      mu1=Voltage/2.0_dp
 	
      write(1,*)Voltage,current(mu1,mu2)  
    enddo
     	    
   end subroutine IvsV_curve
   
!-------------------------------------------------!
! IvsV Single Point 	 	   	 	  !		
!-------------------------------------------------!

   subroutine IVSinglePoint() 
   real(dp) :: G_G0
   integer :: i

   call cpu_time(t_1)     
    
   if(abs(SCPV_muL)+abs(SCPV_muR).NE.0.0_dp)then
     write(*,*)'IV single point ...'    
     singleIV_Current=current(SCPV_muL,SCPV_muR)
     write(*,*)'mu_left=',SCPV_muL
     write(*,*)'mu_Right=',SCPV_muR   	
     write(*,*)'I=',singleIV_Current,'microAmpere'
   endif

   if(abs(SCPV_muL)+abs(SCPV_muR).EQ.0.0_dp)then
     if(Banner_logic.EQ.'S') write(*,*)'0 voltage Fermi level conductance'
     G_G0=transmision(0.0_dp)
     if(WithDecoherence_logic.EQ.'S') call Tdeco1()
     if(Banner_logic.EQ.'S') write(*,*)'G/G0 =', G_G0
     if(Banner_logic.EQ.'S') write(*,*) 'G(nS) =', G0*G_G0
     if(Banner_logic.EQ.'N') write(*,*) G0*G_G0
     if(Banner_logic.EQ.'S') write(*,*)'G Incoherent G/G0', Tdeco
   endif

   call cpu_time(t_2)     
   write(9,*)'IVSinglePoint_time',t_2-t_1
  
   end subroutine IVSinglePoint

!-------------------------------------------------!
! Rho Single Point 	 	   	 	  !		
!-------------------------------------------------!

   subroutine RhoSinglePoint() 
      
   integer :: i
   complex(dp) :: rho_SinglePoint(N,N)
   real(dp) :: total_fill
    
   write(*,*)'rho matrix'
   mu_left=SCPV_muL
   mu_right=SCPV_muR
   call initialize_integration()
   rho_SinglePoint=path_integrator()
   call finalize_integration()
   call write_matrix(rho_SinglePoint)
   
   open(111,file='charges.dat')

   total_fill = 0.0_dp

   do i=1,N
      write(111,*)i,real(rho_SinglePoint(i,i))
      total_fill = total_fill + real(rho_SinglePoint(i,i))
   enddo
      
      write(111,*)'total fill',total_fill

   close(111)
     
   end subroutine RhoSinglePoint

!-------------------------------------------------!
! G(z) Single Point 	 	   	 	  !		
!-------------------------------------------------!

   subroutine GSinglePoint() 
      
   integer :: i
    
    write(*,*)'G matrix'
    mu_left=SCPV_muL
    mu_right=SCPV_muR
    call write_matrix(G_dev(z_SinglePoint))
    
   end subroutine GSinglePoint     

!-------------------------------------------------!
! I vs L courve nanowires 	   	 	  !		
!-------------------------------------------------!

   subroutine IvsL_curve()
      
   integer :: i
    
   do i=4,24
    if(mod(i,2).EQ.0)then
     N=i
!    call h_make()
     call IVSinglePoint()
    endif
   enddo
    
   end subroutine IvsL_curve

!-------------------------------------------------!
! Eigenvals
!-------------------------------------------------!

   subroutine eigenvals()
   
   implicit none
      
   REAL(dp), ALLOCATABLE :: U(:,:),w(:),work(:),Aux(:,:)
   real(dp) :: Homo, Lumo
   INTEGER :: i,info,j,lwork, N_electron
   character :: jobz
   complex(dp), allocatable :: Aux1(:,:)
   
   write(*,*)'calculating eigenvalues ...'
   
   if(UseSMatrix_logic.EQ.'N')then 
     allocate(Solap(N,N))
     Solap=Id()        
   endif 
   
   lwork=8*N-1
    
   allocate(U(N,N))   
   allocate(w(N))
   allocate(work(lwork))
   allocate(Aux(N,N))
   allocate(Aux1(N,N))
    
   Aux=Solap
   Aux1=Solap
   
   U=H	
    
   call dsygv(1,'N','U',N,U,N,Aux,N,w,work,lwork,info)
   
   N_electron = 0

   write(*,*)N_atoms

   do i=1,N_atoms
    N_electron = N_electron + Z_NUMBER(atom_number(i)) 
   enddo

    
   Homo = w(int(N_electron/2.0_dp))
   Lumo = w(int(N_electron/2.0_dp)+1)
   
   if(mod(N_electron-1,2).EQ.0)then  !Paridad
     Homo = w(int(N_electron/2.0_dp)+1) 
     Lumo = w(int(N_electron/2.0_dp)+2) 
   endif

  
   if(calculateFL_logic.EQ.'S') fermi_level = (Lumo + Homo)/2
   
   if(Banner_logic.EQ.'S') write(*,*)'FermiLevel=', fermi_level
   if(Banner_logic.EQ.'S') write(*,*)'LUMO=', Lumo
   if(Banner_logic.EQ.'S') write(*,*)'HOMO=', Homo


!   do i=1,N
!    H(i,i)=H(i,i) - fermi_level   !Esto es lo que hacia antes  que esta mal
!   enddo
   
   H = H - real((fermi_level*Id()).x.Aux1,dp)	   !Esto es lo que hago ahora

   open(333,file='eigenvalues.dat')
   do i=1,N
     write(333,*)w(i) - fermi_level, 10.0_dp   
   enddo
   close(333)   


   if(writeEigenvectors_logic.EQ.'S')then
     open(334,file='eigenvec.bin',form='unformatted',position="rewind")  
     write(334)N  !Esto es necesarion para el programa DOSplot
     do j=1,N
       write(334)U(:,j)
     enddo
     close(334)
   endif


   deallocate(U)   
   deallocate(w)
   deallocate(work)
   deallocate(Aux)
   deallocate(Aux1)


   end subroutine eigenvals
   
!-------------------------------------------------!
! End of module 
!-------------------------------------------------!

   
end
